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ABSTRACT 

Context. The study of physical and chemical properties of massive protostars is critical to better understand the evolutionary sequence which 
leads to the formation of high-mass stars. 

Aims. IRAS 18151-1208 is a nearby massive region (d = 3 kpc, L ~ 2xl0 4 L G ) which splits into three cores: MM1, MM2 and MM3 
(separated by 1'— 2'). We aim at (1) studying the physical and chemical properties of the individual MM1, MM2 and MM3 cores; (2) deriving 
their evolutionary stages; (3) using these results to improve our view of the evolutionary sequence of massive cores. 

Methods. The region was observed in the CS, C 34 S, H 2 CO, HCO + , H 13 CO + , and N 2 H + lines at mm wavelengths with the IRAM 30m and 
Mopra telescopes. We use ID and 2D modeling of the dust continuum to derive the density and temperature distributions, which are then used 
in the RATRAN code to model the lines and constrain the abundances of the observed species. 

Results. All the lines were detected in MM1 and MM2. MM3 shows weaker emission, or even is undetected in HCO + and all isotopic species. 
MM2 is driving a newly discovered CO outflow and hosts a mid-IR-quiet massive protostar. The abundance of CS is significantly larger in 
MM1 than in MM2, but smaller than in a reference massive protostar such as AFGL 2591. In contrast the N 2 H + abundance decreases from 
MM2 to MM1, and is larger than in AFGL 2591. 

Conclusions. Both MM1 and MM2 host an early phase massive protostar, but MM2 (and mid-IR-quiet sources in general) is younger and more 
dominated by the host protostar than MM1 (mid-IR-bright). The MM3 core is probably in a pre-stellar phase. We find that the N 2 H + /C 34 S ratio 
varies systematically with age in the massive protostars for which the data are available. It can be used to identify young massive protostars. 

Key words. ISM: individual objects: IRAS 18151-1208 - ISM: abundances - Stars: formation - Line: profiles 



1. Introduction 



How high-mass stars form is still an open issue (e.g. Zinnecker 



& Yorke 2007). It is particularly not clear whether the for- 
mation process for OB/high-mass stars is different from the 
way solar-type/low-mass stars form. Stars more massive than 
~10 Mq may form like a scaled-up version (high accretion 
rates) of the single (or monolithic) collapse observed for the 
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low-mass stars, or require a more complex process in which 
competitive accretion inside the central regions of a forming 



cluster may play a decisive role (e.g. Bonnell et al. 20041. 
In the first scenario, the observed high-mass clumps (100 to 
1000 M©; 0.5 pc in size) are expected to fragment to form self- 
gravitating high-mass cores (10 to 100 M ; 0.01-0.1 pc in size; 
e.g. Krumholz et al. 2007) which would collapse individually 
to form a massive single or binary stars. In the second scenario, 
the competitive accretion is expected to occur inside the high- 
mass clumps. The study of the properties of high-mass clumps 
is therefore a central observational issue to progress in our un- 
derstanding of the earliest phases of high-mass star formation. 

From a selection of IRAS sources not associated with any 
bright radio source but having the IRAS colors of UC-Hn re- 
gions (as defined by |Wood & Churchwell 1989 1, Sridharan 
|et al.| ( |2002[ ) have built a sample of so-called High Mass 
Protostellar Objects (hereafter HMPOs) which would corre- 
spond to the pre-UC-Hn phase of the formation of high-mass 
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Stars. [BeutheFeFaQ (|2002b} found that the HMPOs were sys- 
tematically associated with massive clumps as detected in the 
dust continuum and in CS line emission. These clumps are 
good candidates to correspond to the earliest phases of high- 
mass star formation. They have not yet formed any bright Hn 
regions and contain a large amount of gas at high densities. The 
precise evolutionary stages of the individual HMPOs might 
however be very diverse and they require to be derived indi- 
vidually through dedicated, detailed studies. 

( 2QQ7\ investigated the whole 



Recently, Motte et al. 



Cygnus X complex and obtained a first unbiased view of the 
evolutionary scheme for high-mass clumps and cores. Like for 
the HMPOs, the massive clumps (~0.5 pc in size) in Cygnus X 
could be resolved into massive cores (~0. 1 pc). Roughly half of 
these cores were found to be very bright in the (mid)-infrared 
(such as AFGL 2591; e.g. |van der Tak et al.|1999| l, and there- 
fore very luminous: the mid-IR high-luminosity massive cores. 
The other half are weak or not detected in the mid-IR, here- 
after the mid-IR-quiet massive cores. Surprisingly all the mid- 
IR-quiet massive cores were however found to drive powerful 
SiO outflows. They could therefore be safely understood as the 
precursors of the infrared high-luminosity massive cores. 

IRAS 18151-1208 is a rath er typical (L ~ 10 4 L Q ) and rel- 
atively nearby (3 kpc) HMPO ( |Sridharan et al.|2002| i. |Beufher 



et al. ( 2002b I show that the clump actually splits into four in 
dividual cores MM1, MM2, MM3 and MM4 (see Fig.[T]for an 
overview of the region). In addition MM1 seems to further split 
into two cores, separated by 16". We hereafter refer to MM1- 
SW for the secondary peak in the South-West of MM1 which 
possibly hosts an embedded protostar ( Davis et al.|20 04). The 
weakest core MM4 is clearly outside the main region, hence 
it will not be further considered in this paper. A CO outflow 
toward MM1 was discovered by |Beuther et al.| ( |2002c| >. The 
MM1 and MM2 cores have sizes roughly two times larger, and 
similar masses (using the same dust emissivity and tempera- 
ture) than the most massive cores in Mot te et al.| l |2007) . The 
IRAS source coincides with MM1 and no significant IRAS 
contribution can be safely attributed to MM2 or MM3. MM3 
is the least massive and compact core and could actually be 
a still quiescent or pre-stellar core. The IRAS 18151-1208 
region is therefore particularly interesting to study since it 
hosts three individual cores which could be interpreted as high- 
mass star formation sites in three different evolutionary stages. 
Millimeter and sub-millimeter wave observations are reported 
in |McCutcheon et al.| ( fT995> and |Beuther et al.| p002b] >. The 
near-IR counterparts (H2 jets and HH objects) of the CO out- 
flow driven by MM1 have been imaged by |Davis et al. ( 2004 1. 

Molecular line observations compared with results of line 
modeling based on a physical model of the source is a classical 
technique to constrain physical and chemical properties of pro- 
tostellar objects (see |Ceccarelli et al.|1996 van der Tak|2005 



for a review). While different approaches can be adopted, the 
most reliable and most often used method consists in constrain- 
ing first the physical (density and temperature distributions) 
model from the dust continuum emission (mid-IR to millime- 
ter wavelengths), and then use this physical model to derive the 
fractional abundances of molecular species from line emission 
modeling ( van der Tak et al.|1999 Hogerheijde & van der Tak 



a [J2000] 

Fig. 1. Image of the MSX 8u.m emission toward 
IRAS 18151-1208 (color scale) overlaid with the 1.2 mm map 
(white polygon) by Beuther et al. (2002b} (white contours 
at 5 and 10 %, and grey contours from 20 % to 90 % of the 
maximum). The large and small crosses indicate the positions 
of the IRAS source and of MM1-SW (see text) respectively. 
The black and white star symbols show the positions of the 
methanol and water masers respectively. The green polygon 
displays the region mapped in the present study (see Fig. |2]i. 



2000; Sc hoier et al.|2002l|Belloche et al.|20021|Hatchell & vail] 
der Tak 2003). Due to the dramatic changes in physical con- 
ditions inside the protostellar envelopes (increase of density 
and temperature, radiation, shocks) chemical evolution is ob- 
served and can be modeled thanks to chemical network codes. 
It is even expected that chemistry could provide a reliable clock 
to date protostell ar objects (e.g. |van Dishoeck & Bla ke 1998; 
Doty et al.|2002j|Wakelam et al.|2004| i~ 



After the description of the molecular observations and of 
the continuum data from the literature (Sect. 2), the results 
are presented in Sect. 3. Section 4 details the modeling proce- 
dure from the fits of the spectral energy distributions (hereafter 
SEDs) for MM1 and MM2 using MC3EQ code (Sect. 4.1.1 
and 4.2.1) to the non-LTE calculations of the line profiles and 
intensities for all observed molecules using RATRAN^ code 
(Sect. 4.1.2, 4.1.3, 4.2.2, and 4.2.3). In Section 5, we discuss 
the results of this detailed analysis in order to improve our ob- 
servational view of the evolution of high-mass cores from pre- 
stellar to high luminosity massive protostars. 



2. Observations. 

Observations were performed during two sessions. The first 
one in June, 2005 at Mopra, the 22m dish of the Australia 



1 See 

2 See 



Wolf et al 



19991 for details. 



van der Tak et al. ( 1999 2000a) for details 
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Table 2. List of observational parameters for IRAM 30m and 22m Mopra telescopes. Here are indicated observed species, 
energy level transitions, line emission frequencies, half power beam width (HPBW), instrument (30m for IRAM-30m and 
Mo for Mopra telescope), main beam efficiency r] m b, receiver name, velocity resolution 6v, observation mode and system 
temperature T svs . 



Species 


Transition 


Frequency 


HPBW 


Instrument'' 


Vmb 


Receiver 


Sv 


Mode 


T 

1 SVS 






[GHz] 


["] 








[rn-s" 1 ] 




[K] 


CS 


7 = 2-1 


97.9809 


25 


30m 


0.78 


A100 


60 


RaM" 


260 




7 = 2-1 


97.9809 


32 


Mo 


0.49 


3.5mm 


191 


OTF r 


240 




7 = 3-2 


146.9690 


17 


30m 


0.69 


D150 


80 


RaM" 


420 




7 = 5-4 


244.9356 


10 


30m 


0.49 


D270 


96 


RaM" 


1150 


c 34 s 


7 = 2-1 


96.4129 


25 


30m 


0.78 


A100 


61 


RaM" 


190 




7 = 2-1 


96.4129 


32 


Mo 


0.49 


3.5mm 


194 


OTF r 


-270 




7 = 3-2 


144.6171 


17 


30m 


0.69 


D150 


81 


RaM" 


390 


CO 


7 = 2-1 


230.5380 


11 


30m 


0.52 


HERA 


406 


RaM2'' 


-600 


H 2 CO 


7 = 3o3 - 2o2 


218.2222 


11 


30m 


0.54 


HERA 


107 


RaM2* 


250 




7 = 3 22 — 2 2 i 


218.4756 


11 


30m 


0.54 


HERA 


107 


RaM2* 


250 




J = 3l2 - 1\\ 


225.6978 


11 


30m 


0.53 


A230 


104 


RaM" 


870 


HCO + 


7=1-0 


89.1885 


36 


Mo 


0.49 


3.5mm 


210 


OTF r 


-260 


H 13 CO + 


7=1-0 


86.7542 


36 


Mo 


0.49 


3.5mm 


216 


OTF r 


-175 


N 2 H + 


7=1-0 


93.1737 


36 


Mo 


0.49 


3.5mm 


201 


OTF r 


-180 



" RaM - Raster Map mode set to create a 3 x 3 pixels mini-map with a A8 step of 10" around source emission peak. 
b RaM2 - Raster Map mode with HERA receiver, set to create a 190" x 140" map with a A9 step of 6" towards MM1, MM2 and MM3. 
c OTF - On-The-Fly observing mode carrying out on a 5' x 5' grid with an angle step Ad of 9" and an OFF position at 30' from the center. 
d Conversional factor is S /T mb = 4.95 ly/K for IRAM 30m telescope and S /T mb = 22 ly/K for ATNF 22m telescope. 



Table 1. IRAS 18151-1208 sources characteristics. Offsets 
from IRAS position, J2000 coordinates and velocity in the lo- 
cal standard of rest ( Beuther et al.|2002b l are reported. 



source 


A«["] 


A6["] 


o-[I2000] 


<5[I2000] 


v [km/s] 


MM1 


13.2 


-4.9 


18'' 17'" 5 8. s 


-12°07'27" 


33.4 


MM2 


-98.9 


-32.8 


18'' 17"' 50.4 s 


-12°07'55" 


29.7 


MM3 


-72.3 


26.5 


18''17"'52.2 ! 


-12°06'56" 


30.7 



Telescope National Facility^] (ATNF). The second one in 
August, 2005 at Pico Veleta, on the IRAM 30m antenna^] 

2.1. Mopra observations 

The IRAS 18151-1208 region was mapped in the On-The- 
Fly (OTF) mode in the rotational transition of CS (7=2-1) 
at 98.0 GHz, C 34 S (7=2-1) at 96.4 GHz , HCO + (7=1-0) at 
89.2 GHz, H 13 CO + (7=1-0) at 86.8 GHz and N 2 H + (7=1-0) 
at 93.2 GHz. Each run covered a field of 5' x 5' with a half- 
power beam width (HPBW) of approximatively 36" (rj m b = 
0.49) and sampled with a AO step of 9". We used the 3 mm re- 
ceiver of the ATNF telescope with the digital auto-correlator 
set to the 64 MHz bandwidth with 1024 channels and both 
polarizations. It provides a velocity resolution of ~0.2 km-s -1 
over a 200 km-s -1 range. Observations happened while weather 
was clear, with a system temperature T sys - 175 - 260 K (see 



Table |2]for details). Pointing and focus adjustments have been 
set on Jupiter or known sources of calibration. 

Raw data have been reduced with AIPS++ LIVED ATA and 



GRIDZILLA tasks ( jMcMulUn et aE1[2004| l. LIVEDATA per- 
forms the subtraction between the SCAN spectra row and the 
OFF spectrum. Then it fits the baseline with a polynomial. We 
used the GRIDZILLA package to resample and build a data 
cube with regular pixel scale, weighted with the system tem- 
perature T sys . We created a 39 x 39 grid of 9" x 9" pixel size 
convolved by a Gaussian with a FWHM of 18", truncated at an 
angular offset of 36". Finally, due to positional errors that we 
found in maps, a shift in a of 14" for N 2 H + , 10" for CS and 5" 
for HCO + has been applied to fit millimeter-continuum peak 



positions of the three sources given by Beuther et al. (2002b). 



3 http://www.narrabri.atnf.csiro.au/mopra/ 

4 http://iram.fr/IRAMES/index.htm 



2.2. IRAM-30m observations 

We observed the rotational transition of CO (7=2-1) molecule 
at 230.5 GHz and the two H2CO rotational para-transitions 
(7=3o3-2o2) and (7=322-22i) at respectively 218.2 GHz and 
218.5 GHz simultaneously, using the HERA 3x3 pixel dual 
multi-beam receiver with VESPA backends at a resolution of 
160 kHz for CO and 40 kHz for H 2 CO. Thus we built two 
190" x 140" maps sampling every 6" around MM1, MM2 
and MM3. This observation was done while weather was good 
{T% m = 0.17), with Tsy S ~ 600 K and 250 K for CO and H 2 CO 
respectively. The maps are incomplete because raster mapping 
was not set to cover the entire Aa and AS, and the HERA ver- 
tical polarization (set for the CO) had a dead pixel (number 4). 
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All pointings and focus have been set on suitable planets (i.e. 
Jupiter) or calibration sources. 

The IRAS 18151-1208 region was also mapped in the rota- 
tional transitions of CS at 98.0, 147.0 and 244.9 GHz (7=2-1, 
7=3-2, 7=5-4), in the isotopic C 34 S rotational transitions at 
96.4 and 144.6 GHz (7=2-1, 7=3-2) and in the H 2 CO rota- 
tional ortho-transition (7=3 12 - 2n) at 225.7 GHz, using si- 
multaneously A 100, A230, D150 and D270 receivers coupled 
to the high-resolution VESPA backend (resolution of 80 kHz). 
Central emission peaks of MM1, MM2 and MM3 were quickly 
mapped with a 3 x 3 raster position switching method with a 
step of 10". Observations happened while sky conditions were 
reasonable (rg' m ~ 0.1-0.6 and T sys ~ 260-1150 K). 

Reduction for all IRAM 30m telescope data were per- 
formed with the CLASS software from the GILD AS suite 
( |Guilloteau & Lucas|2000[ ). We found and eliminated unusable 
data and treated platforming that appeared in CO and H2CO 
baselines, then spectra at the same position were summed and 
finally antenna temperature T* was converted into T„,i, (see 
Table [2] for the values of t] m b). 

3. Results 

The modeling of the spectral energy distribution (SED) of a 
source is a common way to derive its density and temperature 
profiles (e.g. van der Tak et al.| [T999). The obtained so-called 
physical model is then used as such to derive the abundances 
and to further investigate the physical properties (kinematics 
for instance) through the line emission modeling. 

We thus observed CS, C 34 S, HCO + and H 13 CO + line emis- 
sions which are bright and are tracing dense gas. We also ob- 
served H2CO line emission which traces dense gas and its tem- 
perature. In additions we observed N2H + line emission which 
is a cold gas indicator and finally CO that traces molecular gas 
flows. 

3.1. Large scale maps: CO (2-1), H 2 CO (3 22 - 2 2i ), 
CS(2-1), HCO + (1-0), H 13 CO + (1-0) and 
N 2 H + (1-0) 

The velocity integrated maps in CO (2-1), H2CO (3 2 2-22i), 
CS(2-1), HCO + (1-0), H 13 CO + (1-0) and N 2 H + (l-0) are 
displayed in Fig. [2] Table [5] gives the derived emission exten- 
sions (at the 3cr noise level), the minor-to-major axis ratio b/a 
and position angle (RA.) of MM1, MM2 and MM3. Note that 
the C 34 S (2-1) MOPRA map was too noisy and was discarded. 

The molecular emission always peaks at the positions of the 
three cores (millimeter-continuum peak positions by |Beufher| 
et al.|2002b] >. While MM1 and MM2 are always well detected, 
MM3 seems to be too weak in CS (2-1), H 2 CO (3 22 -2 2 i) and 
H 13 CO + (1-0) to be detected. The emission is extended even 
between the cores for CS, HCO + and CO lines but are more 
peaked for N 2 H + , H 13 CO + and H 2 CO lines. Some chemical 
differences are already clear in CS, N2H + and H2CO. While 
MM1 is generally brighter than MM2, such as in the dust con- 
tinuum, MM2 has the same intensity in H2CO and is even sig- 
nificantly brighter in N 2 H + than MM1. 



CS (2-1) line emission is detected for MM1 and MM2 and 
weak emission could be present at l<x for MM3. MM1 emis- 
sion is stronger than MM2, and covers a larger area, with exten- 
sions of 97" x 81" and 75" x 62" for MM1 and MM2 respec- 
tively. It indicates that a large scale envelope is surrounding 
MM1 and MM2. Axis ratios b/a are equal to 0.83 for MM1 
and MM2, showing that the emission is almost spherical. The 
position angle of MM1 emission (+57°) is similar to the posi- 



tion angle of the large dust emission of +62° given by Beuther 



et al. (2002b). We note that emission is detected between MM1 
and MM2. It suggests that they could be connected by a dense 
molecular filament. 

The cold gas tracer N2H + (1-0) is detected toward the three 
cores. Emission from MM2 is stronger and more concentrated 
than for MM1 (89" x 62" and 55" x 54" respectively), sug- 
gesting that MM2 has a smaller and colder envelope. A simi- 
lar behaviour has been observed by |Reid & Matthews] ([2007 ) 
for two similar cores in IRAS 23033+5951. The elongation of 
MM1 (b/a = 0.69) and position angle (+77°) ressembles the 
extension of the dust emission. MM3 line emission is stronger 
compared to other molecules and shows a connection with 
MM2. Together with a similar rest velocity for MM2 and MM3, 
it confirms that MM3 belongs to the IRAS 18151-1208 region. 

HCO + (1-0) and H 13 CO + (1-0) lines are detected for both 
MM1 and MM2. MM3 has a weak emission in HCO + (1-0). 
As for CS, MM1 emission is stronger and more extended 
than for MM2 (respective extensions are 104" x 78" and 
57" x 56"). With b/a = 0.75 and a position angle of +75°, we 
note that MM1 emission may be extended along the outflow as 
expected since the HCO + line emission is definitely showing 
outflow wings such as in CO. The map shows an emission be- 
tween MM1 to MM2, and MM2 to MM3, confirming that the 
three sources are connected. 

Dense gas traced by H2CO (322-2 2 i) line is detected for 
MM1 and MM2 only. The emission from each source coin- 
cides perfectly with the dust continuum positions. MM1 and 
MM2 have similar intensities and have roughly the same small 
circular size (35" x 30", b/a = 0.85 for the two sources) and 
do not show any spatial suggestions of being contaminated by 
outflows. 

The CO (2-1) line is detected in all sources and shows an 
extended emission from the whole region. MM1 emission is 
extended to the west (72" x 51", b/a = 0.71 and PA .= +85°) 
due to the red wing of the outflow (P.A.= +88°). See the fol- 
lowing section for more details on the CO outflows. 

3.2. Molecular outflows in CO (2- 1 ) 

The CO (2-1) line emission shows clear outflow wings. They 
are due to the already known CO outflow driven by MM1 
(Beu ther et al.|[2002c| ) but seem also to be located in the re- 
gion of MM2. Figure [3] displays the average CO (2-1) spectra 
around MM1 and MM2, as well as the contour map of the in- 
tegrated emission in the wings. The displayed spectra clearly 
show emission wings up to high velocities. The observed maxi- 
mum velocities are -20.4 and +20.6 km-s~' for MM1 and -19.7 
and +28.3 km-s -1 for MM2, with respect to vlsr (see Fig. [3] 
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Table 3. Sources extension at 3cr noise , minor-to-major axis ratio b/a and position angle (P.A.; counterclockwise angle relative to 
the North) for each molecular transition observed. When the emission is almost spherical (b/a > 0.85) P. A. is uncertain and set 
to0°. 







MM1 






MM2 






MM3 




Ext. 


b/a 


P.A. n 


Ext. 


b/a 


P.A. n 


Ext. 


b/a 


P.A. ["] 


CO (2-1) 


72" x 51" 


0.71 


+85 


34" x 23" 


0.68 


-53 


45" x 36" 


0.8 


+53 


H 2 CO(3 22 - 2 21 ) 


35" x 30" 


0.85 





35" x 30" 


0.85 











CS (2-1) 


97" x 81" 


0.83 


+57 


75" x 62" 


0.83 


-31 








N 2 H + (1-0) 


89" x 62" 


0.69 


+77 


55" x 54" 


0.98 





54" x 52" 


0.96 





HCO + (1-0) 


104" x 78" 


0.75 


+75 


57" x 56" 


0.98 











H 13 CO + (1-0) 


52" x 38" 


0.73 


-28 


< HPBW 













■ 12°06'0(J" 



-ia°06 r oo" 



■12°06'00" - 



-12°0B'00" 



-12"09 , 00" 




lfi h lfl m 00 s 



Fig. 2. Integrated velocity maps of CS (2-1), N 2 H + (1-0), HCO + (1-0), H 13 CO + (1-0) obtained with the Mopra 22m telescope 
(left and middle panels) and CO (2-1), H2CO (322 _ 22i) maps obtained with the IRAM 30m telescope (right panel). Level con- 
tours fit 90 % to 10 % of peak flux with a 10 % step. Only levels exceeding the threshold of 2cr noise in each map are plotted in order 
to identify detections. Peak fluxes are 16.9 K km s 1 for CS, 16.9 K-km-s" 1 for N 2 H\ 32.9 Kkms 1 for HCO + , 2.31 Kkm-s' 1 
for H 13 CO + , 21.1 K-km-s -1 for H2CO and 393 K-km-s -1 for CO. Crosses indicate the 1.2 mm continuum positions (Beuther et al. 
2002b) for the three sources. 



and Table HJ. It is clear from the map of this wing emission 
that there are actually two outflows driven by the two respec- 
tive cores. The MM2 outflow is actually newly discovered (see 
Sect. 5.1 for a discussion). 

The MM1 outflow appears mostly east- west oriented. This 
is not quite the orientation observed for the main jet-outflow 
system recently discussed by Davis et al. (2004) with a NW-SE 
direction for the H2 jet (see their Fig. 2). A second H2 jet, west 
of MM1 and probably driven by an embedded source in MM1- 
SW (see Fig. [TJ, is however also recognised in this work. The 



position of MM1-SW is marked in Fig.|3]with a small green tri- 
angle. Its location is close to a possible secondary center of out- 
flow with a weak blue-shifted lobe in the south and the bright 
red-shifted western lobe being curved toward that location. The 
spatial resolution of our CO observation is not high enough to 
fully conclude about the association of CO outflowing gas with 
the different driving sources, but it seems that the western part 
of the MM1 outflow might actually be dominated by a second 
outflow from MM1-SW. The inclination of the MM1 flow on 
the sky is difficult to determine. From the significant extension 
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of the flow and from the overlap of the lobes we can only say 
that the flow is neither face-on (in the plane of the sky) nor 
pole-on, and is therefore mostly intermediate. 

The MM2 outflow is more compact than for MM1. Both 
lobes are confined in the core. The integrated emissions in the 
two lobes are similar in strength and are four times larger than 
the blue lobe of MM1. If the flow is dominated by a single 
ejection, the large overlap of the two lobes suggests that the 
flow is not face-on, but could be close to pole-on. 

All the derived properties of the two outflows are summa- 
rized in Table |4] In addition, to quantify the energetics of the 
two detected outflows, we derive and give in this table the out- 



flow forces using the procedure described in Bontemps et al. 
( 1996 ). The mechanical force (or momentum flux) seems to be 
the measurable quantity from observed CO outflows which can 
approach in a not too uncertain way the corresponding quan- 
tity of the inner jet or wind ( jCabrit & Bertout|[T992} |Chernin| 
|& Masson|1995| l. We estimate the outflow force Fco for MM1 
and MM2 by integrating the momentum flux inside volumes 
with radius 6" and 18" centered on the sources. In addition, 
for the extended red lobe of MM1, a measurement is given be- 
tween 18" and 30" to better cover the red lobe. Since the incli- 
nation of the two outflows can not be reliably derived from such 
low resolution observations we adopt the average correction of 



a factor of 10 adopted by Bontemps et al. ( 1996). 



3.3. Raster maps in CS, C 34 S, and H 2 CO (3 n -2 n ) 

In addition to the large maps, a complementary investigation 
of the three millimeter sources is provided by the 3x3 raster 
maps of CS, C 34 S and H 2 CO (see Fig.|4j). 

All transitions are detected toward MM1 at u^sr = 33.4 + 
0.1 krn-s . The line profiles are similar with a typical line 
width (FHWM) Au = 2.7 ± 0.2 km-s" 1 . CS(2-1), CS (3-2) 
and CS (5-4) have roughly equal intensities. The optically thin 
C 34 S(2-1) line emission is detected all around the source at 
an approximately constant intensity level of ~ 1 K. In contrast 
C 34 S (3-2) seems to peak on the central position. The line pro- 
files and the spatial patterns for H2CO (3i2-2n) and CS (5-4) 
are very similar, suggestive of common origin and excitation 
conditions for the two transitions (E\xp/k are respectively equal 
to 33.40 and 35.27 K). 

The intensities of most transitions toward MM2 are weaker 
than MM1. The line profiles at u LSR = 29.8 ± 0.1 km-s -1 have a 
triangular shape and a larger line width {Au = 3.7 ±0.1 km-s -1 ) 
compared to MM1. This could be due to a gas motion contri- 
bution through its turbulence. The intensity of CS (5-4) is low 
compared to MM1, contrary to H2CO (3i2-2n) line emission 
which is as bright as in MM1. 

All molecular line profiles are weaker in MM3, and seem 
to be less peaked than in MM1 and MM2 as seen in velocity 
integrated maps (cf. Fig.[2|. They are detected at all positions in 
CS (2-1), (3-2), they have a weak emission in two positions in 
H 2 CO(3i2-2n) and finally CS (5-4) and the C 34 S transitions 
could not be detected. 
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Fig. 3. (a) Contour maps of the CO (2-1) blue-shifted (dashed) 
and red-shifted (solid) integrated emission in the outflow wings 
for MMI, MM2, and MM3 (no detection). Positions of the 
MM cores are marked by filled (green) triangles. Blue- and 
red-shifted velocity intervals are +13.0 to +27.4 and +39.4 
to +54.0 km-s- 1 for MMI, +10.0 to +23.8 and +35.8 to 
+58.0 km-s -1 for MM2 and MM3. Lowest, highest contours 
and contour steps are (10, 70, 15) K-km-s -1 (blue) and (60, 410, 
50) K-km-s- 1 (red) for MMI, and (10, 70, 15) K-km-s-' (blue) 
and (20, 195, 25) K-km-s-' (red) for MM2 and MM3. The 
background grey-scaled image is the integrated H2CO emis- 
sion from Fig. ||] (b) CO (2-1) spectra for the MMI outflow 
together with the velocity cuts used to integrate maps (dotted 
vertical lines) in (a) and to derive the mechanical forces. The 
solid line circles around MMI and MM2 are 6"and 12" ra- 
dius and indicate the area on which the Fqo measurements are 
derived (see text for details). The filled grey spectrum is an av- 
erage over the region of bright intensities in the wings as seen 
in the integrated map in (a). The dashed and blue (respectively 
solid and red) spectrum shows the line profile at the maximum 
of the blue (resp. red) lobe, (c) same as (b) for MM2. 
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Table 4. Momentum flux of CO (2-1) calculated as in Bontemps et al. 
blue-shifted momentum flux between 6" and 12" Fco b (M Q -km-s^'-yr 



19961. Blue and red integration intervals Aii£ >r (km-s 



(M Q -km-s 
MM2. 



), red-shifted momentum flux between 6" and 12" Fqo, 
•yr _1 ), total momentum flux Fco lMlll (M -km-s^'-yr -1 ), maximal and outflows extensions are reported for MM1 and 



Table 5. Continuum flux densities of the three cores within IRAS 18151-1208 region. MSX flux densities are extracted from 
the catalog for MM1 and a 3>cr noise calculation over the beam is done to derive an upper limit for the undetected MM2 and MM3 
sources. IRAS measurements are derived from maps improved by HiRes method. Millimeter and sub-millimeter continuum flux 
densities are adapted to source reference size derived from continuum map by Beuther et al. (2002b). 



MSX flux densities [Jy] IRAS flux densities [ly] mm-continuum flux densities [Jy] mm source 
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33.0 


61.8 


19.0 


98.6 


<891 


<1890 


49.4 
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MM2 
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<0.34 


<0.73 


<0.38 


<1.36 
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0.670 


18.1 


14.4 


MM3 


<0.07 


<0.26 


<0.41 


<0.82 


<0.61 


<2.21 


<891 


<1890 








0.261 


18.0 


14.2 



3.4. Continuum fluxes 

In order to derive the spectral energy distribution of the three 
millimeter peaks in IRAS 18151-1208, we retrieve continuum 
data from the literature. MM1 was imaged for the first time 
at 450, 800 and 1 100 ] xm by |McCutcheon et al] ( [1995) then 
at 450 and 850 \xm by | Williams et al.| (|2004]i7 and the over- 
all region at 1.2 mm by |Beuther et al. (2002b). Uncertainties 



are homogenized to account for typical calibration errors taken 
as 20 % in the millimeter range (1.1 and 1 .2 mm) and 30 % in 
the submillimeter range (450 and 800 |im). In order to compare 
observations and model, we choose the deconvolved dust emis- 
sion size 6 r as a reference for the SED. We deconvolve using 

(P\ where Q s is the observed mm-continuum half- 



peak size derived by Beuther etaL| ( p00 2b ) and mb the HPBW 
(11"), using a two-dimensionnal Gaussian fit, assuming that 
the source is smaller than the beam and has a Gaussian shape. 



All millimeter (A > 450 ii.m) peaks S , 
this reference size using 



peak 



are then adapted to 



Concerning mid- and far-infrared flux density measure- 
ments, MSX detected MM1 only. Thus we derive a 3<x noise 
level inside the beam size around the millimeter-continuum po- 
sition of MM2 and MM3, taken as an upper limit. We note that 
the IRAS emission of the region is dominated by MM1 which 
is confusing upper limits for MM2 and MM3. We therefore im- 
proved the IRAS map using the publicly available HiRes maxi- 
mum correlation method ( Au mann et al.|1990] >, showing MM1 
detection at 12 \xm and 25 \im. At 60 \im and 100 \im, a sin- 
gle peak is detected slightly shifted to the west. We conclude 
that the huge beam width of IRAS sums the three flux densities 
from the three sources and the environmental filament contri- 
bution. Measured flux densities reported in the IRAS catalog 
at 60 [im and 100 \xm are taken as upper limits for the three 
sources. All results are shown in Table|5] 



4. Modeling the continuum and the molecular 
emission of MM1 and MM2 



(1) 



et al 



where p is the density power-law index (derived by Beuther 
2002b I and S x the adapted flux density used to build 



the SED. These flux densities are 49.4, 7.62, 1.40 and 1.66 Jy 
at respectively 450, 850, 1100, 1250 |a.m for MM1. The pro- 
cess is iterated on MM2 and MM3 and reference sizes with 
millimeter-continuum peak fluxes are reported in Table [5] 
Unfortunately the sub-millimeter continuum data for MM2 and 
MM3 are less abundant than for MM1 as these two sources 



were not observed by |McCutcheon et al. ( 1995 ) and Williams 
|etaL| ( |2004) . 



In this section we explain how we derive physical models 
for MM1 and MM2 and how we use them to investigate the 
kinematics and the molecular abundances inside the cores. 



Following the method initiated and validated by van der Tak 
|et al.| ( |1999| |2000b| l we use the dust continuum emission to 
derive the density and temperature profiles of the cores. Dust 
emission is sensitive to the total column density of dust (op- 
tically thin part of the SED in the millimeter wave range), 
and to the dust temperature and its distribution(mid to far-IR 
range). We try to build the simplest model in order to have the 
most general representation of the object. Then from the ob- 
tained physical model, we can run molecular line modeling to 
constrain kinematics (line profiles) and chemistry (abundances 
of observed species). Since MM3 is weak or not detected in 
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Fig. 4. 3 x 3 maps in CS, C 34 S, and H 2 CO (3i 2 -2„) for MM1, 
MM2 and MM3 sources. Velocity ranges from 15 to 50 km-s -1 
in all spectra, intensity ranges from -2 to 5 K for H2CO line 
emission, -1 to 2 K for C 34 S line emissions, -2 to 5 K for 
CS (5-4) line emission and -2 to 10 K for other CS line emis- 
sion. Spectra have been smoothed to an identical velocity reso- 
lution of 0.2 km-s -1 . 



4.1. IRAS 18151-1208 MM1 
4.1 .1 . Spectral energy distribution 

The MM1 core is massive (M « 550 M according to 
Beu ther et al.||2002b"| and applying the correction of 1/2 from 
the associated 2005 erratum). It coincides with the bright 
IRAS 18151-1208 source (L * 2.0 x 10 4 L ). In first approx- 
imation, assuming that the object has a spherical symmetry, 
we want to derive the best radial distribution of density and 
temperature that can fit the SED and that takes into account 
observational constrains. 

(1) ID modeling 

We use the radiative transfer code MC3D in its ID version 
(Wolf et al.| 1 999) > to calculate the expected dust emission. This 
code is Monte-Carlo based, and is using the standard MRN 
( Mathis et al.|1977|l dust grain di stribution with the dust opaci- 
ties from Ossenkopf & Henning ( 1994). The model parameters 



are the central stellar luminosity L t , its temperature T t , the in- 
ternal and external radii ro and r ext of the object, the molecular 
hydrogen density no at ro and the power-law index p of the 
density distribution of the form n(r) = rto(r/ro) p . 

The bolometric luminosity is derived by integrating the 
SED, and we obtain L/,,,/ = 1.4 x 10 4 L . The temperature 
of the central star is taken as the corresponding main sequence 
star, i.e. T* = 22 500 K (type B2V of 10 M approximatively, 
cf. Underhill et al. 1979| l. We adjust ro to fit the dust subli- 
mation radius where T ^ 1500 K. The external radius r ext of 
the model is derived from the radius of the 1.2 mm contin- 
uum emission (here the deconvolved mean of major and minus 
axe given by Beuther et al.||2002b i assuming that the source 
has a Gaussian shape smaller than the beam size. The result- 
ing extension is 18.2", hence r„, = 27 300 AU or ^ 0.13 pc 
at 3 kpc. The power-law index p is set to -1.2 according to 
the continuum profile fit by Beuther et al. (20025}. Finally the 



only free parameter left is no and is derived by fitting optically 
thin 1 .2 mm continuum emission with a x 2 calculation. We 
find n = 2.8 x 10 9 cirT 3 at r = 21 AU with X 2 = 3.94, 
hence (n) = 8.0 x 10 s cirT 3 and a total mass of gas equal 
to 660 M Q . Outer and mean temperatures, formally T ext and 
(T), are respectively 25.4 K and 27.3 K. The derived tem- 
perature profile can be fitted with a power-law of the form 
log 10 (D = a. log r + p with a = -0.588 and /? = 3.95. 

The resulting SED (see Fig. [5] in dashed line) shows 
that millimeter and sub-millimeter continuum parts are well 
modeled. On the contrary, the mid- and far-infrared emission 
is not reproduced, due to spherical symmetry that does not 
permit infrared photons, coming from inner parts, to leave the 
source without being reprocessed into photons at millimeter 
wavelengths. 2D modeling, where the polar regions contain 
less matter than in the mid-plane, should therefore reproduce 
better the observed mid-infrared emission. 



a number of molecules and does not have a well constrained 
SED, we do not attempt to model it. We will restrict ourself to 
a simple discussion of its evolutionary stage. 



(2) 2D modeling 

We use the 2D version of the MC3D code which was de- 
veloped for disks. It is parametrized in radial and altitude di- 
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Fig. 5. Spectral energy distributions obtained from ID model 
(dashed line) and 2D model (continuous line) overlaid on ob- 
served fluxes of MM1 source. Dust continuum peak fluxes 
have been adjusted to fit radial extension of the model. Mid- 
infrared fluxes are obtained in 2D choosing a typical view angle 
~ 45° (90° = edge-on view). 

rections (r,z) and it assumes the Keplerian disk equilibrium by 



Shakura & Syunyaev ( 1973 1: 



n(r, z) 



(2) 







otherwise 



We note that the model has three free parameters: «o deter- 
mined as for the ID modeling, the inverse height of scale c, and 
the source angle relative to the line of sight 8i„ s . It is mostly 0/ OJ 
which controls the fraction of mid-infrared emission. Values 
close to 45° can roughly reproduce the SED, and we do not 
further adjust the value to keep the model as general as pos- 
sible. The inverse height of scale c is controlling how thin the 
disk-like structure is. Since MM1 is young and is still mostly a 
spherical core, we try to keep c as small as possible while hav- 
ing always a significant mid-infrared emission for typical aver- 
age lines of sight. Thus c was gradually increased from 0.4 by 
0.1 steps until an infrared emission greater than 0.1 Jy (limit of 
detection in the observations) could be obtained, finally leading 
to c — 0.7. Other parameters are kept identical as in ID model 
(T*, L m , ro, r ext and p). Calculations show that the best fit is 
obtained for «o = 2.7 x 10 9 cirT 3 with^- 2 = 0.732 (total mass 
M = 830 M Q ). The derived SED is plotted in Fig.[5]in continu- 
ous line. It illustrates how a 2D modeling can reproduce better 
the observed infrared emission, although the silicate feature is 
not matched. We display the resulting density and temperature 
distributions in Fig. [6] One can see that the equatorial regions 
are denser and colder than around the polar axis. This is why 
mid-IR emission can escape and directly contribute to the SED. 

4.1 .2. CS and C 34 S line emissions 

Since molecular line emission is usually dominated by the 
mostly spherical external layers of the protostellar cores, the 




r/r„ [%] 



r/r„ [%] 



Fig. 6. Density and temperature distributions for the 2D model 
presented in a logarithmic grey scale, with 20 logarithmic lev- 
els from 10 2 to 10 7 cirT 3 for density and 10 logarithmic levels 
from 10 K to 263 K. 



2D description might not be required to reproduce molecular 
lines. 

We use the CS and C 34 S line emission to compare the re- 
sults of modeling in ID and 2D geometries. We directly use 
outputs of continuum models described above, assuming that 
dust and gas temperatures are equal, to which we add param- 
eters for molecular abundances and kinematics. The molecu- 
lar abundances relative to H2, X(CS) and X(C 34 S), mostly de- 
termine the line intensities when the opacity is low (r < 1). 
Then we use a constant turbulent velocity vj- It is the main 
contributor to the line width compared to the thermal disper- 
sion. X(CS) and vj are therefore the free parameters to fit the 
intensities and widths of the line emission we observed. Finally 
we want to test if an infall motion in our model could im- 
prove our results, giving at least an upper limit for it. The in- 
fall velocity is parameterized b y the free-f all, power-law dis- 
tribution v in (r) = v in (r/r Q y - 5 (|Shu||l977|l. We test its effect 



on modeled line emission by switching between u,-„ = and 
v in = -0.4 km-s -1 , the lowest value that modifies modeled line 
emission significantly for MM1. 

The line profiles are computed with the RATRAN code 
(Ho gerheijde & van der Tak||2000| l. This code allows ID and 
2D modeling using respectively multi-shell and grid descrip- 
tion of the object. The radiative transfer calculation uses the 
Monte-Carlo method to derive the populations of the energy 
levels and then builds data cubes that are convolved with corre- 
sponding beam sizes of the telescope to be directly compared 
with the observations. 

In the ID (spherical) geometry, a x 1 calculation is per- 
formed for the optically thin C 34 S line (t 2 -i = 9.0 x 10~ 2 and 
T3-2 = 2.2 x 10- 1 ) over a grid for / C 34 S = X(C 34 S)/X(C 34 S) 07? 
(with X(C 34 S) typ = 1.0 x 10- 10 ) and v T . The best fit is obtained 
for (/ C 34 S ,ur) = (0.5, 1.0) with X 2 = 0.956. Then the grid is 
refined for / C 34 S only, because model results are less dependent 
of ut (kept equal to 1.0 km-s -1 ). The best result is obtained 
for /c ,4 s = 0.5 1 with x 2 = 0.952, leading to the abundance 
X(C 34 S) = 5.1 x IP" 11 . Assuming a typical [S/ 34 S] isotopic 
ratio of 20 dWilson & RoodffT994l |Chin et al.|[T996l |Lucas &| 
|Liszt|[l998] l the CS emission is modeled with the abundance 
X(CS) = 1.0 x 10~ 9 . 

The final results of the line modeling are displayed in Fig. [7] 
Given the very small number (2) of free parameters, the result- 
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ing fits are quite good. It reproduces well the line intensities. 
Only the profiles for the higher excitation lines indicate a sig- 
nificant difference between observations and this simple model. 
There are clear self-absorptions for the (3-2) and the (5-4) 
transitions (rs_4 = 3.6) which indicate an overestimate of large 
CS column densities for the highest excitation regions (central 
regions). This could well be due to a significant decrease of the 
CS abundance toward to the center of the core (see Sect. |5.5| for 
a discussion). 

We find that a minimum infall velocity of -0.4 km-s -1 is 
needed to start to see a clear asymmetric shape of the pro- 
file. Introducing infall improves only slightly the results for the 
(2-1) and (3-2) transitions and makes the fit of the (5-4) tran- 
sition worse (cf. Fig. [7J second plot from the left). We conclude 
that an infall velocity in the ID model does not improve the fits, 
and we consider as an upper limit the value of -0.4 km-s -1 . 

For the 2D modeling, the same^- 2 adjustment is performed 
for vj and /o^s- Interestingly enough, the final best fit is exactly 
the same than in the ID geometry with X(CS) = 1.0 x 10~ 9 
and X(C 34 S) = 5.1 x 10". We note that fine profiles and 
relative intensities between the transitions are similar to the ID 
case. Only a very slight improvement of the profile of CS (5-4) 
which is less self-absorbed can be noted. The result is shown 
on the third plot from the left in Fig [7] 

Again the effect of an infall velocity on the line profiles is 
tested. With i>,„ = -0.4 km-s -1 we note that the asymmetric 
profile of the line emission is less strong than in the ID case 
(cf. Fig.[7J fourth plot from the left). We however conclude that 
infall in the 2D case does not improve critically the results. 

Comparison between ID and 2D modeling shows that us- 
ing a spherical description of the source is sufficient to re- 
produce CS and C 34 S molecular line emission, even if 2D is 
slightly better with less self-absorption than in the ID case. In 
addition we conclude that infall motion is not clearly detected 
and is not necessary to reproduce observations of MM1. We 
thus decide to use only ID modeling without infall motion for 
the other observed lines. 

The fact that 2D modeling is not needed for the CS line 
emission can be explained by the optical depth which is not 
so high at the corresponding frequencies compared with the 
mid-IR continuum, or by the asymmetry which is significant 
only on small scales where the temperature is high (between 
300 and 1000 K). The observed molecular lines are not tracing 
these inner regions. 

4.1.3. Other lines: N 2 H + , HCO+, H 13 CO + and H 2 CO 

We adopt the ID model derived in the previous section. For 
each molecule, the only free parameter left is therefore its 
relative abundance to H2. 

(1) N 2 H + (7=1-0) modeling 

The N2H 4 " (7=1-0) is split by hyperfine structure. 
Theoretically there are 15 hyperfine components, which blend 
into seven for sources with low turbulence (Cas elli et al.|1 995) 
or into three if turbulence is strong, as in our case. In order 
to correctly fit observations from ATNF telescope at Mopra 



we created the molecular datafile for the RATRAN radiative 
transfer code. Energy levels, Einstein coefficients A„/ and colli- 
sional rates y u \ are derived from Daniel et al. ( 2004 ) through the 
BASECOL database maintained by the Observatoire de Paris^] 
Energy levels included in the molecular datafile vary from 
7 = to 7 = 6. The hyperfine structure resulting from angular 
momentum and nuclear spin interaction (7"i and F quantum 
numbers) leads to a total of 55 energy levels. The statistical 
weight of each level is determined by (see Daniel et al.|2005) : 



n J,Fi,F 



(27 + 1)(2F + 1) 
[J,F U F] 



(3) 



where [7, F\ , F] is the number of magnetic sub-levels for the 
angular momentum 7. The molecular data file takes into ac- 
count 119 radiative transitions where the first 15 of them fit to 
the (7=1-0) transition. The spectra obtained are then summed 
to finally reproduce the whole composite profile of the triplet 
(see Fig. [8]). 

Collision rates y u \ reported by Daniel et al. ( 2004} are given 
for helium as collision partner. As described by Schoier et al.| 



(2005) we take 1.4 x y u \ in order to correct for an H2 colli- 
sion partner. Initially collision rates are given for temperatures 
from 5 K to 50 K which is the temperature range of the outer 
parts of our sources. To be able to model the entire sources, we 
have extrapolated the rates to high temperatures using a least- 
squares method that derives collisional rates from a linear fit of 
the form: 



log(r»/) = a.log(T) + b 



(4) 



This extrapolation has the advantage of fitting the global trend 
of the collisional rate without any exaggerated values for high 
temperatures. 

Emission of N2I-L is considered as optically thin (t ~ 10~ 2 ) 
due to its typical low abundance. Thus we model it by com- 
paring the velocity integrated fluxes of observation and model 
starting from a typical abundance X(N2H + ) = 1.0 x 10~". 
The ratio of modeled to observed fluxes is then used to further 
correct the abundance while verifying that the line profile 
is always correctly reproduced. Iterating the process gives 
X(N 2 H + ) = 3.5 x lO" 10 . An overlay of the model on the 
observation is shown in Fig. [8] 

(2) HCO + and H I3 CO + (7=1-0) modeling 
The HCO + and H 13 CO + lines are good tracers of high density 
gas and of its kinematics. H 13 CO + (7=1-0) is usually optically 
thin (here we finally get r = 3.4 x 10~ 2 ) so we adopt the 
same routine to derive molecular abundances as for N2H + . 
We then get an abundance of X(H 13 CO + ) = 3.4 x 10" 11 . 



Applying a typical isotopic ratio [ 12 C]/[ 13 C] = 67 ( jWilson &| 

line 



Rood 1994 



Lucas & Liszt 1998]) we model the HCO + 

'. Results are shown in Fig. [9] 
line is optically thick (r = 1 .7) 
4 K. It suggests that the 



with X(HCO + ) = 2.3 x 10- 
The resulting modeled HCO + 
reaching a maximum intensity of 
region of emission in the model is smaller than the telescope 
beam as expected for a high density tracer. On the other 
hand the observed intensity and profile are very different. It 

http://amdpo.obspm.fr/basecol/ 
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Fig. 7. Emission of CS and C 34 S from MM1 object (continuous line) overlaid by -from left to right- ID static, ID with infall 
(v in - -0.4 km-s -1 ), 2D static and 2D with infall (i>,„ = -0.4 km-s -1 ) model line emissions (dashed). Spectra velocity resolutions 
are kept identical to observational resolutions reported in Table [2] 
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Fig. 8. Emission of N 2 H + from MM1 object overlaid by ID 
model of the source (dashed line). The triplet shape is com- 
pletely reproduced using our RATRAN molecular datafile. 
Spectra velocity resolutions are kept identical to observational 
resolutions reported in Table [2] 



speaking, HCO + is a good outflow tracer and we speculate that 
the HCO + line is actually dominated by some HCO + rich gas 
associated with the outflow shocks inside the core (i.e. at high 
density as required to excite HCO + ). This would mean that 
the derived HCO + abundance in the core envelope is only an 
upper limit and could well be much smaller, except locally in 
outflow shocks. 

(3) para- and ortho-H2CO modeling 
The ratio between para and ortho populations of H2CO is an in- 
teresting parameter in the global context of gas temperature and 
chemical activity (Kah ane et al.|19 84) and can be a useful tool 
to follow chemical evolution of gas inside and between the ob- 
served sources. Our modeling indicates that the lines are mostly 
optically thin (t 2 i8.2ghz = 1-37, t 2 i 8 .5GHz = 0.09 and t 22 5GHz = 
0.77) so we adopt the same routine as described for N 2 H + . 
The abundances obtained are as follows X(para - H 2 CO) = 
1.5 x 10 -10 and X(ortho - H 2 CO) = 1.3 X 10 -10 , leading to a 
ratio para/ortho ^1.2. The agreement between the modeled and 
the observed lines is good. All the lines however tend to show 
some excess on the blue side, while the ortho line even shows 
excesses on both sides (see Fig. 10 1. In addition we want to 



emphasize that the two para-H 2 CO transitions which have two 
significantly different Eup/k (respectively 21.0 K and 68.1 K) 
are well reproduced. This is an indication that the excitation 
temperature of the emitting gas is well reproduced by the phys- 
ical model. 



indicates that the bulk of the HCO + rich gas does not follow 
the general distribution of matter indicated by other tracers 
(dust emission and optically thin lines). The HCO + (7=1-0) 
line seems to have an excess on the blue side which could 
well be associated with the blue-shifted outflow wing which 
is mostly located inside the core in CO (see Fig.|3]l. Generally 



4.2. IRAS 18151-1208 MM2 
4.2.1 . Spectral energy distribution 

Since MM2 is not detected in the mid-infrared, and taking into 
account our conclusions from the modeling of MM1, we decide 
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Fig. 9. Emission of HCO + (top) andH 13 CO + (bottom) (7=1-0) 
from MM1 source overlaid by ID model of the source (dashed 
line). Resulting abundances are X(HCO + ) = 2.3 x 10~ 9 and 
X(H 13 CO + ) = 3.4 x 10~ n , according to the typical isotopic 
ratio [ 12 C]/[ 13 C] = 67. Spectra velocity resolutions are kept 
identical to observational resolutions reported in Table [2] 
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Fig. 10. para-H 2 CO (top and middle) and orfho-H 2 CO (bot- 
tom) emissions from MM1 source overlaid by ID model 
of the source (dashed lines). Resulting abundances are 
X(para-H 2 CO) = 1.5 x 10~ 10 and X(ortho - H 2 CO) = 
1.3 x 10~ 10 , leading to a ratio para/ortho =s 1.2. Spectra ve- 
locity resolutions are kept identical to observational resolutions 
reported in Table [2] 



to only model the source in a simple ID, spherical geometry. 
We cannot derive the bolometric luminosity of MM2 by inte- 
grating its SED, due to the fact that all infrared fluxes are upper 
limits. Instead, we assume a luminosity of 2700 L , which is 
uncertain by at least a factor of 2. This luminosity corresponds 
to a B4V type star (approximatively 7 M G ) with T„ = 16 600 K, 
that we use to describe the heating source of our model. The 
1 .2 mm continuum deconvolved source extension from Be ufherl 
et~aT] ( (2002bl > is equal to 14.4" thus r ext =21600 A.U. or 
= 0.10 pc at 3 kpc and we assume a power-law index p = —1.3 
from the same paper. Because only one mm-wave peak flux 



density was measured for MM2 (see Table [5]), we adapt no 
to fit this unique measurement. After two iterations we con- 
verge to no = 1.1 x 10 10 cirT 3 at r = 12.4 A.U., hence 
(n) = 1.1 X 10 6 cm~ 3 and a total mass gas of M - 460 M . We 
derive T ext = 19.2 K, (T) = 19.4 K and the temperature can be 
fitted by a power-law of the form log 10 (r) = a. log r + f3 with 
a = -0.614 and /3 = 3.85. The resulting SED does not show 
any emission in the mid- and far-infrared, in agreement with 



the observations (see Fig. 1 1 1 and as expected for a simple ID 
description. 
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Fig. 11. Spectral energy distribution of MM2 overlaid with the 
ID model result. Infrared emission has always an upper limit 
due to the non-detection of the source. 

4.2.2. Molecular line emission modeling 

We use the same procedure to model the molecular lines to- 
ward MM2 as for MM1. We therefore have derived the molec- 
ular abundance relative to H2, and the turbulent velocity vj. 
We have also tested a possible infall velocity field of the form 
Vin(r) = v in (r/r )-°- 5 . 

Emission of optically thin C 34 S lines (t 2 _i = 4.4 x 10~ 2 and 
r 3 _ 2 = 4.4 x 1(T 2 ) in MM2 object are treated by taking two dif- 
ferent grids for the x 2 calculation over its intensity and area. 
The first grid is (/ C 34 S ,i;r) = ([0.1, 0.2... 0.5], [1.5, 1.6.. .2.0]) 
and shows a best fit for (f,v T ) = (0.3, 1.7), the second one 
is finer with / C 34 S = [0.25, 0.26. ..0.35] and a fixed vj - 
1.7 km-s -1 because of the low dependence of the modeled line 
on this parameter, compared to fc^s- It gives the best fit for 
/ C 34 S = 0.27, hence an abundance X(C 34 S) = 2.7 x 10" 11 . As 
for MM1, the resulting modeled lines reproduce well all tran- 
sitions except the CS (5-4) line which is heavily overestimated 
in intensity. Again it points to a lower abundance of CS in the 
inner regions compared to the outside. 

Then, as above for MM1, we add an infall velocity profile 
following a typical r~ 05 distribution. The asymmetric shape is 
obtained for a minimum infall velocity i>,„ = -1.0 km-s -1 . The 
resulting profile is unchanged for (2-1) transition, is slightly 
improved for the (3-2) but shows an excess of blue emis- 
sion which is not observed for CS (5-4) (see Fig. 12 right). 



Therefore including infall velocity distribution does not im- 
prove our fit and hence does not seem to be necessary. It gives 
an upper limit for the infall velocity of u,„ = -1.0 km-s -1 . 

The optically thin emission of N 2 H + (r ~ 10~ 2 ) in MM2 
is treated by comparing the observed and modeled velocity in- 
tegrated area and using it to adapt an initial typical abundance 
(1.0 x 10" 10 ). The resulting fit for N 2 H + is shown in Fig. [l3 
and corresponds to X(N 2 H + ) = 6.3 x 10~ 10 . As for MM1, the 
quality of the fit is high. 

For HCO + and H 13 CO + we obtain the following abun- 
dances: X(HCO + ) = 5.1 x 10" and X(H 13 CO + ) = 7.6 x 10" 11 . 
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Fig. 12. Emission of CS from MM2 object overlaid by ID 
static (left) and ID with infall (y,„ = -1.0 km-s -1 ) model of 
the source (dashed lines). Spectra velocity resolutions are kept 
identical to observational resolutions reported in Table [2] 
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Fig. 13. Line emission of N 2 H + from MM2 object overlaid by 
the ID model of the source (dashed line). Spectra velocity res- 
olutions are kept identical to the observational resolutions re- 
ported in Table [2] 



vations. The HCO + column density derived from the H 13 CO + 
line implies an heavily saturated (r = 19) HCO + line which is 
not observed. Obviously in contrast to all the other molecules, 
the distribution of HCO + does not follow the simple spherical 
geometry of the model. 

We then derive the best abundances for the para and ortho- 



The best fit result is shown in Fig. 14 As for MM1, the mod- 
eled HCO + line emission does not reproduce well the obser- 



H 2 CO to reproduce the three observed transitions (see Fig. 15 1. 
The abundance obtained for the ortho transition at 225 GHz 
is Z(ortho - H 2 CO) = 2.0 x 10" 10 . For the two para tran- 
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Fig. 14. Emissions of HCO + (top) and H 13 CO + (bottom) (7=1- 
0) from MM2 source overlaid by ID model of the source 
(dashed line). Spectra velocity resolutions are kept identical to 
observational resolutions reported in Table [2] 
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Fig. 15. para-H2CO (top and center) and ortho-H2CO (bot- 
tom) emissions from MM2 source overlaid by ID model of 
the source (dashed lines). Spectra velocity resolutions are kept 
identical to observational resolutions reported in Table [2] 



sitions, no unique abundance can be derived. For the best fit 
of the 218.222 GHz line emission, the abundance obtained is 
X(para-H 2 CO) = 2.6 x 10" 10 but the 218.475 GHz line is 
then too weak by a factor of almost 2; see Fig. [15] At the other 
extreme, if the 218.475 GHz line is fitted, the abundance ob- 
tained is X(para - H 2 CO) = 4.8 x 10~ 10 . In contrast to MM1 
for which a correct ratio of the two para lines was obtained, the 
modeled ratio (218.222 over 218.475) for MM2 appears to be 
too large. Since the weaker 218.475 GHz line emission has a 
higher upper energy level (Eup/k = 68.1 K) it suggests that 
the emitting gas is actually warmer in MM2 than what is rep- 
resented by the simple ID model. The resulting para to ortho 
ratio is therefore in the range =s 1.3- 2.4. 



5. Discussion 

5.1. MM2: a new massive protostar driving a powerful 
outflow. 

MM2 is a millimeter source without infrared counterpart (MSX 
or IRAS). It is even seen in absorption over the local back- 
ground in the MSX 8iim image (see Fig. [I}. A wat er maser 
has been detected toward the core by |Beuther et al. (2002a I. 
With the addition of our newly discovered CO outflow driven 
by MM2, it all clearly points to a protostellar nature of MM2. 



Inside the core size of 0.22 pc (14.4", see Table 3.4 1 a total 
mass of 460 M is obtained with MC3D (dust emissivity for 
the MRN grain distribution; see Table [6j. MM2 is therefore 
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" Values from continuum map analysis by Beuther et al. (2002b). 

Table 6. Summary of results from dust continuum emission 
modeling (above mid-line) and from molecular line emission 
modeling (below mid-line). The relative uncertainties on the 
derived parameters are given inside the parentheses and are 
at the 3 cr level. They correspond to the rms noise propagated 
through the modeling process. The additional absolute uncer- 
tainty on masses and abundances, mostly due to the uncertain 
emissivity of dust, is given in Sect. |5.3| 



a newly discovered massive protostellar core which is not de- 
tected in the mid- or far-IR: an IR-quiet massive protostellar 
object. A more complete comparison with the detailed analysis 



obtained for low-mass protostars and in Cygnus X by Motte 
et al. ( 2007) is given in the following section. 



5.2. Evolutionary stages of MM1, MM2 and MM3. 



In contrast to low-mass star formation (e.g. An dre et al.|20 00), 
the evolutionary sequence for high-mass stars is not well con- 
strained and understood. The general lack of spatial resolution 
leads to discussions of evolutionary sequences mostly applied 
to massive clumps (such as HMPOs, typical size of 0.5 pc) or 
to high-mass dense cores (typical size of 0. 1 pc) which cannot 
usually be directly compared to protostars which have physical 
sizes more of the order of 0.01 - 0.05 pc (e.g. |Motte et al.|1998| 



Mott e & Andre|200T i. A recent attempt to clarify these differ- 
ent scales and associated evolutionary sequences is compiled in 



Motte et al. (2007 1; see their Table 4. From the complete sur- 
vey of the relatively nearby Cygnus X complex, a total of 40 
high-mass dense cores were recognized with an average size of 
0.13 pc, similar to sizes of nearby, low-mass dense cores but 
20 times more massive and 5 times denser. These high-density 
cores were all found to be already protostellar in nature and 



were proposed to well represent the earliest phases of high- 
mass star formation. Among the 40 cores, 15 were already as- 
sociated with an UC-Hn region and bright in the infrared, 8 
were only bright in the infrared, and 17 were IR-quiet proto- 
stellar cores. This sequence is proposed to be mostly an evo- 
lutionary sequence, the IR-quiet cores being the precursors of 
the IR-bright cores, and then of UC-Hn regions. We will dis- 
cuss the cores in IRAS 18151-1208 in the framework of these 
recent Cygnus X results. 

First of all no bright radio source is detected in the whole 
IRAS 18151-1208 clump. All the massive protostellar objects 
detected in the clump are necessarily objects in an evolutionary 
stage prior to the formation of an UC-Hn region. In Fig. [T] the 
distribution of mid-infrared emission from MSX is displayed. 
It is clear that the main bright infrared source is associated with 
MM1 while the MM2 area is devoid of mid-infrared emission. 
The core is actually even seen in absorption over the local back- 
ground. MM3 is more complicated since a moderately bright 
MSX source is situated close to the center of the core. 

MM1 has a mass of 660 M G in a radius of 27 300 A.U. (core 
size of 0.27 pc) obtained from the fit of the SED with MC3D 
(Table It has a bolometric luminosity higher than 10 4 L , 
and is bright in the mid-IR (S20 nm = 61 Jy). It drives a powerful 
bipolar outflow and therefore hosts at least one mid-IR bright 
massive protostar. Using the same dust properties (millime- 



ter dust emissivity and average temperature) as in Motte et al. 
(2007), and inside a 0.13 pc size, the MM1 core would have a 
mass as high as 74 M G , and would therefore be among the most 
massive cores of Cygnus X. With an equivalent S20 ^ m = 190 Jy 
at 1 .7 kpc (distance of Cygnus X), it would be the brightest core 
which is not associated with an UC-Hn region (see Fig. 7 in 
Motte et al.|2007 1. It is however not as bright in the mid-IR as 



the well-known AFGL 2591 source described by |van der Tak| 
|eT1u^ ( fT99"9"] l. In the |Motte et^L] ( p007] > classification, MM1 is 
therefore a high-luminosity IR (or IR bright) protostar. 

Except its non detection in the mid- and far-IR (S20 Hm 5 
2 Jy), the properties of MM2 ressemble those of MM1 . It has a 
mass of 460 M in a radius of 21 600 A.U. (core size of 0.22 pc) 
obtained with MC3D (Table|6|, a still high bolometric luminos- 
ity of 2700 L Q , and it drives a powerful bipolar outflow. Using 
the same dust properties as in Motte et al7| ( |2Q07| >, and inside a 
0. 13 pc, the MM2 core would have a mass of 45 M Q , and would 
therefore be among the 40 high-density cores of Cygnus X. In 
the |Motte et aL] ( |2007[ ) classification, MM2 is a mid-IR-quiet 
massive protostar. 

Observations of MM3 does not permit to conclude defini- 
tively on its nature. Although it is massive with a flatter density 
profile (about 200 M Q and p = 0.8, tjBeufher et al.|2002b) and 
no outflow detected, any mid-IR emission is confused due to 
the nearing confusing source detected by MSX and the pres- 
ence of a bright IR filament in the background. At best one can 
only propose that MM3 is a probable prestellar core. 

One can also use the CO outflow energetics to further dis- 
cuss the evolutionary stages of MM1 and MM2. Outflows are 
believed to be the best indirect tracers of protostellar accre- 
tion. Interestingly enough, while MM2 is not a bright IR source 
and is less luminous than MM1, its CO outflow is as power- 
ful as the one of MM1. This is like Class low-mass proto- 
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stars which have on average more powerful outflows than more 
evolved, and IR bright Class I objects. It has been interpreted 
by Bontemps et al. ( 1996| > as due to a decrease of the accretion 
rate with time. Extrapolating these results to higher masses, we 



investigate in Fig. 16 how MM1 and MM2 place in the evolu- 
tionary diagram based on the energetics of CO outflows. The 
location of the low-mass protostars (Class and Class I) is dis- 
played with stars and is well reproduced by a toy model (dot- 
ted tracks) based on an exponential decrease of the accretion 



rate with time (see details in Bontemps et al.| [1996) for stars 
between 0.2 and 2 M . Three additional tracks are given for 
M = 8, 20, 50 M Q . One can note, for instance, that while 
the low-mass protostars reach a maximum luminosity which 
is higher than their final ZAMS luminosity (by a factor of 5 
for 0.2 M Q stars) the luminosities of massive protostars are al- 
ways increasing. In this scenario, the dashed line would repre- 
sent the location of the transition between Class and Class I 
protostars, i.e. the location where half of the final mass is ac- 
creted (first black arrow symbols, second arrow symbols are for 
90 % of the mass accreted). If this can be extrapolated to high- 
mass protostars it indicates that MM2 is younger than MM1 
and could host a high-mass Class like protostar. Note that 
the HMPO (star symbols) including MM1 and MM2 are not 
well resolved as individual collapsing objects. This is partic- 
ularly important for MM2 for which the measured luminosity 
depends on its continuum fluxes in the millimeter range and 
therefore directly scales with the size of the region. The lumi- 
nosity of MM2 should therefore be seen as an upper limit of 
the actual luminosity of the protostar. In this diagram, MM2 
and MM1 are located in the area of the M = 20 M track with 
MM2 having less than 50 % of its final mass, and MM1 hav- 
ing of the order of 50 % of its final mass accreted. While the 
stellar embryo of MM2 is therefore certainly not yet a ionising 
star, MM1 could be hosting a weakly ionizing stellar embryo 
of ~10 M Q . This could explain why MM1 would appear as 
the brightest core in Cygnus X (in Fig. 7 of Motte et al.|2 007 ) 
which is not yet associated with a UC-Hh region. 

The overall properties of MM1 and MM2 therefore provide 
us with a coherent picture of MM1 and MM2 hosting at least 
one high-mass protostar which are not yet creating an UC-Hn 
region. MM2 hosts a certainly younger protostar, and proba- 
bly slightly less massive, than MM1. The high-mass protostar 
hosted by MM1 is a good candidate to be just at the limit of 
switching on a new UC-Hn region. In this scheme, well-known 
protostars such as AFGL 2591 would be slightly more evolved 
with an already formed UC-Hn region and a much brighter flux 
in the mid-IR. 



5.3. Density and temperature profiles from dust 
continuum 

From the detailed ID and 2D analysis of the continuum emis- 
sion and of the optically thin molecular line modeling, we ar- 
rived at the conclusion that a ID spherical geometry is enough 
to reproduce the molecular line emission. It is only to account 
for the mid-infrared emission that a 2D geometry is required. 
Our careful analysis showed that the ID and 2D models were 
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Fig. 16. Diagram of the bolometric luminosities versus CO 
outflow momentum fluxes for protostars adapted from Fig. 5 
in Bontemps et al. ( 1996 ». Big starred markers (black and 



white): low-mass protostars (Bontemps et al. T996JI; small 
starred markers: unresolved massive protostars in HMPOs 
from Beuther et al. (2002c). The dotted tracks arrows are from 
the [Bontemps et al.| ([1996} toy model for masses from 0.2 to 
50 M Q . The black circle for MM1 indicates the values obtained 



by Beuther et al. ( 2002c i for this source, the black squares the 



values we derived and given in Table |4]and[6] The white square 
indicates the location of the MM2 high-mass protostar. The to- 
tal luminosity of the MM2 core is taken as an upper limit of the 
luminosity of the unresolved high-mass protostar. 

indistinguishable when it comes to model the molecular lines. 
This is actually not surprising since only a very small fraction 
of the total mass in the core is affected by the 2D inner flatten- 
ing of the density distribution. Only 1 1 M Q of the 660 M of 
MM1 is found to be at T > 100 K. The bulk of mass is at 
much lower temperatures and does not contribute to the in- 



frared emission (see Fig. 17 I. 

Facing the same difficulty to reproduce the mid-infrared 
emission of massive protostars, |van der Tak et aT] ([1999) 
adopted another strategy by using a ID spherical model with a 
large inner cavity which could let escape the infrared photons. 
Our results also validate this approach since it was a way to ne- 
glect the inner dust contribution for the molecular line model- 
ing. On the other hand, this approach tends to converge toward 
large sizes of the cores which are not observed. We believe it 
is more realistic to impose the sizes of the cores as observed in 
the dust continuum. 

Finally, we note that the use of a self-consistent dust radia- 
tive transfer code such as MC3D is ideal to reduce the num- 
ber of free parameters. On the other hand, it imposes some 
stringent hypotheses on the dust properties by requiring a full 
opacity curve for dust grains, from infrared to millimeter wave- 
lengths. For instance, we use the MRN distribution of grains 
which is known not to take into account very well the real dust 
opacity in the millimeter range when ice mantles are present 



for cold and dense cores ( |Ossenkopf & HennTn g 1994). For 
SED modeling, this is not very important since dust emission 
in the millimeter range is optically thin. But it is however im- 
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Fig. 17. Density (top) and temperature (bottom) profiles used 
in ID for molecular line emission modeling of MM1 (plain) 
and MM2 (dashed). Intervals (vertical lines) represent where 
90 % of the mass is present in the models hence what physical 
conditions are dominating line emission process of molecules 
observed. 



5.4. Average abundances from molecular 
observations. 

Modeling of MM1 and MM2 allowed us to derive absolute 
molecular abundances. The CS abundances we find, X 0.5- 

1.0 x 10~ 9 , are typical values obtained for this type of source 
( Pirogov et al.|2007 1 . We note the same r esult for N2H + where 
X ^ 3.5-6. 3 x lQ-' u (|Pirogovetal.|2003t , HCO + with X * 2.3- 

5.1 x 10~ 9 ([van der Tak et al.|2000a), and finally X(H 2 CO) is in 



the large typical abundance range observed ( |Keane et al.|2001| 
Bisschop et al. 20071. All these results fit abundance ranges 
de rived from stand ard chemical evolution modeling as made 
by |Doty et al] ( |2002| : X(CS) varying from lxlO" 10 to lxl0~ 8 , 
X(N 2 H + ) from lxlO" 12 to 3xl0~ 10 , X(HCO + ) from lxlO 11 to 
6xl0~ 9 and X(H 2 CO) from lxlO" 10 to 2xl0~ 8 . Thus we con- 
clude that our observations and models do not reveal any abun- 
dance anomaly compared to other massive protostellar objects 
and to chemical predictions. 

5.5. Depleted CS abundance in the inner regions. 

The modeling of multiple line emission from a single molec- 
ular species enables us to probe physical conditions inside a 
source. In our model, where no shocks are inserted, the main 
factor is thermal excitation. Modeling of CS transitions (£\ip Ik 
is equal to respectively 7.0, 14.1 and 35.3 K for 7=2-1, 7=3-2 
and 7=5-4) reveals that for both MM1 and MM2 sources the 
two lowest CS transitions are almost reproduced, whereas the 
highest one is stronger than observed (see Fig.[7jand Fig. 12 1. A 
lower abundance of CS in the inner regions could explain this 
discrepancy. Since the gas is cold and dense in basically the 
whole MM1 and MM2 cores a significant depletion of CS onto 
the grains is expected (see Tafalla et al. 2002, and references 
therein). For the line modeling of the low-mass Class proto- 
star IRAM 04191, |Beflbche et al.| ( [2002] ) had also to consider a 
significant CS depletion (by a factor of about 20) to reproduce 
the observed lines. The CS depletion is expected to be still sig- 
nificant even for high-mass cores such as MM1 and MM2 be- 
cause they have not yet formed a large hot core region (only a 
very small fraction of the total mass of MM1 is expected to be 
at T > 100 K; see above). We conclude that CS depletion in the 
inner regions of MM1 and MM2 could be responsible for the 
observed too weak CS(5-4) emission which is not reproduced 
by our modeling. 



portant for mass derivation, and we therefore expect that the 
masses we derived are overestimated due to freeze out onto 
grains. As a consequence the abundances could be underesti- 
mated accordingly. The resulting absolute uncertainty can be 
evaluated by directly deriving total masses from the 1 .25 mm 
continuum fluxes of MM1 and MM2 in Table fusing the mass 
determination by Motte et al. ( 2001) and the average temper- 
atures obtained in Table [6] We get 200 and 114 M for MM1 
and MM2 respectively which are ~3.5 times smaller than with 
MC3D. This systematic effect is however a general concern and 
our resulting abundances are still directly comparable with re- 
sults from most previous similar studies which used the same 
assumptions. 



5.6. H 2 CO: a probe of physical and chemical 
interactions ? 

Our study clearly shows that H 2 CO emission of we observed 
can be modeled with no radial variation of abundance. The test 
of an overabundance driven by a hot core (T > 100 K) does 
not show any significant improvement in line fitting. We can 
understand it from the low contribution in mass of the hot parts 
of MM1 and MM2 (resp. 1 1 M G and 3 M G ; see Fig. [17]), Thus 



we conclude that H 2 CO line emission does not show any hot 
core contribution. 

However the line emission from MM2 leads to a greater 
abundance of H 2 CO than MM1, with a need for a hotter inner 
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part although it is the colder and the younger of the two sources. 
The line emission from MM2 could be affected by a significant 
contribution of outflow shocks in the inner envelope of the pro- 
tostar. This idea is supported by the higher turbulent velocity 
observed in MM2 (vj = 1.7 km-s" 1 in contrast to 1.0 km-s -1 
for MM1) and its internal para-to-ortho ratio greater than 1, 



suggesting recent chemical activity ( Kahane et al. 1984 1 driven 
in these shocks. 



5. 7. Chemical evolution of high-mass cores. 

We finally wish to speculate that the chemical differences be- 
tween MM1 and MM2 are related to their different evolution- 
ary stages. The derived abundances, obtained with the same 
modeling process and radiative transfer code, show that X(CS) 
is 2 times higher for MM1 than for MM2, and 5 times lower 
than in the probably more evolved APGL 2591 high-mass pro- 
tostar, suggesting that the CS abundance might increase with 
time. Interestingly enough, the Johnstone et al. ( 2003| > study 
of sub-millimeter sources in Orion was suggesting the same 
trend. The chemical models by Wake lam et al.| ( (2004) > actually 
predict such an evolution of the CS abundance in protostars: 
a constant production of CS is expected at low temperatures 
(T<100 K) thus increasing the abundance with the source evo- 
lution. Moreover, when the hot core region develops, the re- 
lease of CS from the grains may even increase more the CS 
abundance. This could be the case for the high abundances in 



AFGL2591 and in the sources observed by Hatchell & van der 
rMl ( |2003] l. 

In contrast, X(N2H + ) is found to decrease with time; the 
N2H + abundance inside AFGL 2591 shows a clear trend of 
decrease from MM2 to AFGL 2591: X(N 2 H + ) AF gl 2591 = 
0.1 x X(N 2 H + )mmi = 0.05 x X(N 2 H + ) M M2- This molecule 
hardly depletes on cold dust grains and is rapidly destroyed 
at warm temperatures. This behaviour is not unique and has 



been already observed in high-mass protostellar objects (Reid 



& Matthews 2007). We suggest that the CS and N 2 H + , i.e. 



abundance ratio may be a good tracer of protostellar evolution, 
but more observations and modeling are required to test this 
hypothesis. 

6. Conclusions 

Here we summarize our conclusions on the massive protostellar 
objects MM1 and MM2 of the IRAS 18151-1208 region. 

1. The three cores of the region, MM1, MM2 and MM3 are 
physically linked and have probably been formed from a 
single parental cloud or clump. 

2. We detected for the first time a CO outflow driven by MM2. 
It clearly establishes the protostellar nature of MM2. In 
contrast MM3 does not show any outflow activity and is 
therefore most probably a pre-stellar core. 

3. Following the evolutionary scheme discussed in Motte et al. 
(2007), MM1 is a high luminosity IR (or IR-bright) massive 
protostar while MM2 is an IR-quiet massive protostar. 

4. We have established that while an inner flattening of the 
matter distribution is required to reproduce the SED of 



MM1, a simple ID spherical geometry is enough to well 
model molecular line observations. In contrast, MM2 does 
not even require a inner flattening since it is not detected in 
the mid-infrared range. 

5. A significant depletion of CS in the inner parts of the MM1 
and MM2 cores is required to fully reproduce the observed 
CS line emission. 

6. We find that the abundance ratio between CS and N 2 H + 
could be a very good evolution tracer for high-mass proto- 
stellar cores hosting high-mass protostars. 
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